Screened empirical bond-order potentials for Si-C 
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Typical empirical bond-order potentials are short ranged and give ductile instead of brittle behav- 
ior for materials such as crystalline silicon or diamond. Screening functions can be used to increase 
the range of these potentials. We outline a general procedure to combine screening functions with 
bond-order potentials that does not require to refit any of the potential's properties. We use this 
approach to modify Tersoff's [Phys. Rev. B 39, 5566 (1989)], Erhart & Albe's [Phys. Rev. B 
71, 35211 (2005)] and Kumagai et al.'s [Comp. Mater. Sci. 39, 457 (2007)] Si, C and Si-C poten- 
tials. The resulting potential formulations correctly reproduce brittle materials response, and give 
an improved description of amorphous phases. 

PACS numbers: 



I. INTRODUCTION 

Empirical and semi-empirical classical interatomic po- 
tentials have been used in computer simulations for more 
than two decades. Bond-order potentials (BOPs) 1 - 
a class of semi-empirical formulations — have proven 
to yield reasonably accurate potential energy landscapes 
for covalently bonded 2-19 and metallic 20-30 materials. 
The bond-order approach can be systematically derived 
from the tight-binding approximation. 8,20 ' 21 ' 31 ' 32 This 
furnishes the hope that although simple, BOPs should 
show transferability to a wide number of situations. 

The rigorous derivation of BOPs by Pettifor and co- 
workers 8 ' 31,32 was predated by empirical formulations 
that are the scope of this article. 2 ' 3 ' 33,34 The parameters 
of these empirical BOPs are adjusted to match ground- 
state properties, such as the cohesive energies or elastic 
constants, for one or more phases of an element or com- 
pound. For covalently bonded materials, the interaction 
between atoms is usually limited to nearest-neighbors 
and also limited to short distances. Both limitations are 
independent of each other, although a short interaction 
range is typically used to limit the interaction to nearest- 
neighbors. This is possible because in crystalline struc- 
tures second-nearest neighbors are well-separated from 
first neighbors. They show up as well-distinguishable 
peaks in the atomic pair distribution functions. It is less 
clear that limiting the interaction range works in liquids 
or amorphous solids where such separation is not nec- 
essarily given. Additionally, the short range makes the 
description of transition events such as the dissociation 
of a bond inaccurate. Qualitatively unphysical behavior 
is obtained in particular when a transition is driven by 
external forces. 1,35,36 

An example of this latter problem is a crack that is 
driven through a brittle material. In contrast to phys- 
ical reality, empirical BOPs consistently predict duc- 



tile behavior for materials such as silicon or carbon. 35 
This problem is usually circumvented by using full quan- 
tum calculations 3 ,38 or by embedding a quantum region 
around the crack tip in a classical potential. 39 ' 40 How- 
ever, the interaction between multiple cracks, driving a 
crack in an amorphous material, or a series of mode II 
cracks such as a tribological interface 41-47 would be noto- 
riously difficult to model with either approach. A classi- 
cal interatomic potential that reproduces brittle fracture 
is therefor highly desirable. 

We have shown in an earlier work 36 that brittle 
behavior can be restored by decoupling the condi- 
tion for nearest-neighbor relationship from the range 
of the potential. This potential was based on the 
second-generation reactive empirical bond-order poten- 
tial (REB02), 10 and nearest-neighbor relationship was 
determined using the screening functions first introduced 
by Baskes and co-workers 48 in the context of the mod- 
ified embedded atom method. This modification kept 
the REB02's ground-state properties of crystalline struc- 
tures and molecules untouched. Here, we add two silicon- 
carbide potentials and one pure silicon potential to the 
family of screened BOPs. The first two are based on Ter- 
soff's 3rd 4 and Erhart & Albe's 13 potentials. The Si-only 
potential is based on the parameterization by Kumagai 
et al. that was optimized for the melting point of sili- 
con. 49 The particular form of Kumagai's potential has a 
desirable feature that Tersoff's and Erhart & Albe's are 
lacking. All potentials are modified in a manner that does 
not change the properties of crystalline ground states. 



II. SECOND-MOMENT BOND-ORDER 
POTENTIALS 

In bond-order potentials of the Tersoff-Brenner type, 
the cohesive energy E of a structure is expressed as a sum 
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over bonds. For each bond the energy has a purely re- 
pulsive, <p(r), and a purely attractive, /3(r), contribution. 
The strength of the attractive contribution is modulated 
by the bond-order, a quantity that depends on the en- 
vironment of the bond and that is related to Coulson's 
bond-order concept. 1 ' 50 The particular expression we use 
here is 



E=lj2S ij [<P(r ij )-b ij /3(r ij )] 



(1) 



i<j 



where 4>{rij) and j3(rij) are pairwise positive functions. 
The function bij is the bond-order. SV, is a switching 
function that switches the interaction off under certain 
conditions and will be described in more detail below. 
The expression for the bond-order is 



Refs. 4 and 13 uses trigonometric functions and is given 

by 



fc(r) = < 



1 + COS ( 7T 



\ r 2 -ri J 



if r < r\ 

if n < r < r 2 (8) 
if r > r 2 



The interatomic potential that is defined by Eqs. (1) to 
(8) with minor differences in the choices of 4>{r) and j3(r) 
and bij has been used to parameterize, among others, the 
interaction of B-C-N, 7 ' 9 C-H, 5 - 6 ' 10 C-O, 12 C-O-H, 18 Ga- 
As, 11 Fe-C, 27 < 29 Pt-C, 22 Si-C, 2 - 4 - 13 Si-C-H, 19 Si-O, 16 - 17 
W-C-H 24 and Zn-O. 25 While it is possible to go beyond 
second moments to higher chemical accuracy, such po- 
tentials have only been developed for few element com- 
binations, such as Mo 23 W 28 Fe 30 C-H 8 < 14 and Si. 15 



with 



*« = (1 + X&)- 
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III. SMOOTHNESS OF INTERATOMIC 
POTENTIALS 



Xij 



E 



Sikh(rij,r ik )g(Oij k ) 



where 



h(rij,r ik ) = exp{[2/z ife (r y - r ik )] m } . 



(3) 



(4) 



Here, g{9) is some function with angular periodicity and 
rj, S, fi and m are free parameters. The bond-order en- 
ables a directional dependence of bonding and hence sta- 
bilizes the open cage-like structures (e.g. diamond) that 
covalently bonded materials form. 

In what follows, we will discuss potentials where the 
pairwise functions <j> an d /? are given by exponentials of 
the form 



/3(r) = exp 



and 



0(r) 



Dn 



K-l 



exp 
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-av2K(r — r ) 



(5) 
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With these choices, the dimer potential energy curve 
(3(r) — <j>(r) has its minimum at ro with energy Dq and 
curvature a. The angular function is typically given by 



, c 2 c 2 
9(0) = 7 ( 1 + ^ - d 2 + [h + cos6} 2 



(7) 



where c, d and h are free parameters. The parameter K 
characterizes the relationship between equilibrium bond 
energy and bond length for different crystal structures. 

The range of all distance-dependent functions is lim- 
ited to nearest-neighbors by a switching function Sij = 
fc( r ij ) that depends on the distance between atoms i and 
j only. Generally, fc drops from a value of one to zero 
between two distances r\ and respectively. A com- 
mon choice of cut-off function that is also employed in 



Experience tells that potential energy surfaces as ob- 
tained for example from density functional theory (DFT) 
calculations are smooth. This fact is underlined by the 
recent success of using Gaussian processes 51 to extrap- 
olate from a finite set of energies obtained from DFT 
calculations to arbitrary configurations. 52 A couple of 
DFT calculations typically suffice to reconstruct high- 
accuracy potential energy surfaces. In Gaussian pro- 
cesses, smoothness is intrinsically programmed into the 
extrapolation by the covariance function. 

The potential energy landscape obtained from Eqs. (1)- 
(8) is not smooth because the cut-off function Eq. (8) 
forces energies to zero within a short distance interval. 
This leads to a failure in the description of transition 
states that is most easily demonstrated for the dimer. 
Fig. 1 shows the energy and tensile force of the carbon 
dimer as computed using Tersoff's and Erhart & Albe's 
potential. The energy drops to zero steeply as the cut- 
off is approached. This leads to an overestimation of the 
force required to break this bond, with implications for 
the simulation of cracks and tribology. We here general- 
ize the meaning of the switching function Sij with Sij = 1 
meaning that a bond exists. This allows to unlock the 
asymptotic behavior that is programmed into <j>(r) and 
P(r) for any structure. 

As already mentioned above, the cut-off function is 
designed to allow interaction of nearest neighbors only. 
Physically, this can be motivated by the fact that the 
bond-integral (here the attractive part /? of the poten- 
tial) follows a different functional form for second and 
farther neighbors that is smaller in magnitude. 53 The in- 
teraction of second and farther neighbors is "screened" 
by the nearest-neighbor atom. A central approximation 
in empirical bond-order potentials is to assume perfect 
screening for second and farther neighbors and set their 
bond integral to zero. In a tight-binding (molecular or- 
bital) picture, the physics of screening functions can be 
traced back to non-orthogonality. 54 
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FIG. 1: Energy E and tension / = dE/dr as a function of dis- 
tance r for the carbon dimer. Shown are results obtained us- 
ing the original formulation of the potentials and the screened 
version (denoted by +S) that unlock the asymptotic behavior 
of the interaction. 



A. Cut-off procedure 

Besides finding nearest neighbors, a cut-off criterion 
needs to be able to smoothly interpolate upon transi- 
tions that involve changes in coordination number. We 
have recently proposed to determine nearest-neighbor re- 
lationship 36 from the screening function introduced by 
Baskes et al. 48 that fulfills this condition. Later, Kuma- 
gai et al. 55 have proposed an almost identical scheme. 

The procedure is a follows: Instead of counting atoms 
within a certain distance towards a bond, we look for 
third atoms in the vicinity of the bond. If any third 
atom sits close to the bond it is screened, if it sits far 
away, the bond is allowed to persist. In this picture a 
bond is unscreened if there is a line of sight between the 
two atoms participating in the bond. A simple empir- 
ical and quantitative measure for this intuitive picture 
is given by constructing ellipsoids of revolution through 
two atoms. If a third atom sits inside these ellipsoid the 
bond is screened. 

Let rij denote the distance between atom i and atom j 
for which we would like to compute whether interaction 
is possible. We construct an ellipsis through a third atom 
k (see Fig. 2a). With Xn- = (rik/rij) 2 the coefficient 



2(Xik + Xjk) — (Xik — Xjk) 2 — 1 
1 - (X l k - X jk ) 2 



(9) 



gives the square of the ratio of the two half axes' lengths. 
We now consider a bond between atoms i and j to be 
entirely screened by atom k if the coefficient falls below 
a critical value C m ; n , while an unscreened bond corre- 
sponds to Cijk > C max . A geometric explanation for the 
Cijk coefficient is given in Fig. 2. 





FIG. 2: (a) Screening of the bond i-j: If an atom k moves into 
the vicinity of the bond i-j we construct an ellipsis through 
atoms i, j and k such that bond i-j constitutes one of the half- 
axes. The square root of the coefficient Cijk is then the ratio 
of the lengths of second half axis to first half axis. The Cijk is 
an empirical measure for how close atom fc sits to bond i-j. (b) 
Our screening approach distinguishes two cutoff radii. The 
bond i-j always exists if < n. For r-z < we compute 
the screening function of panel (a) and determine whether 
bond i-j exists from the value of Cijk- If an atom sits in the 
inner gray region where r\ < rij < we interpolate between 
the screened and unscreened cutoff function (see Eq. (11)). 
For reasons of computational efficiency we furthermore turn 
the interaction completely off for r\j > r\ ■ 



We now impose the cutoff on the value of Cijk rather 
th an rij. We define the screening function Sjj of bond i-j 
to be given by Sjj = if the bond i-j is entirely screened 
and otherwise by 48 



n 

k , C j jh^ Crr 



exp 



ijk. 



a 



ijk 



(10) 



The product runs over all atoms k which are neighbors 
to the bond i-j. For each neighbor k we test whether 
atom k might screen the bond, and multiply the contri- 
butions to the screening function accordingly. Addition- 
ally, we do not want the screening to be active in high 
pressure situations, where solids may be compressed to 
highly coordinated structures. Hence, we define an inner 
core region where screening is inactive by choosing the 
switching function to be (see Fig. 2b) 



Sij = fsinj) + (l - f s {nj))T,ij 



(11) 



Here fs is a function that drops from unity to zero be- 
tween radii r\ and r 2 where we switch from a bond that 
cannot be screened to a bond that can be screened by 
its neighbors. Note that is differentiable more than 
twice. To make the overall potential energy landscape 
differentiable more than twice we use: 



fs(r) = { 



1 

exp 




\ r 2 —ri j 



if r < r\ 

if n < r < r 2 (12) 
if r > r 2 



This switching procedure does not introduce an addi- 
tional (artificial) length scale and is intrinsically infinitely 



4 



ranged. The "infinite range" is manifested by the fact 
that all distances occurring in Eq. (9) are normalized by 
the bond distance r«. 



B. Long-ranged limits of the bond-order term 




The long-rangedness necessitates an additional mod- 
ification to traditional empirical bond-order potentials. 
The switching function S appears in the total energy 
Eq. (1), but also in the definition of the bond-order 
Eq. (3). Since for most potentials fi = we find 
h(rij,rik) — 1 and hence the bond-order bij becomes in- 
dependent of the actual bond length and approaches the 
wrong limit in some situations. One of these situations 
occurs when a crystal is cleaved to expose two surfaces. 
As we pull the crystal apart to introduce two free surfaces 
the total energy of the system needs to asymptotically 
approach the energy of two separated systems. 

For the specific bond i-k shown in Fig. 3 the bond- 
length increases continuously with increasing separa- 
tion x. The values of 4>(ri k ) and /3(ri k ) then drop to zero 
as rik — > oo. However, bond i-j feels the presence of 
atom k in the three body term b^. For h = 1, this term 
is given by 



bij = 




(13) 



and independent of the absolute length ri K of bond i-K if 
that bond is unscreened and Si K = 1. For the particular 
bond i-k shown in Fig. 3 we have Si k = 1. Without 
any mechanism to eliminate the influence of atom n = k 
to the bond-order in bij in Eq. (13) the bottom surface 
will feel the top surface's presence at arbitrary distances 
since Sik = 1- Without screening functions we have Sik = 
fc( r ik) depend only on distance, and the contribution of 
k will have vanished once the atom has moved out of the 
cut-off radius of atom i, i.e. once n k > Ti. 

Here we argue that in order to provide a well defined 
limiting value for the bond-order with increasing bond- 
length we must choose fi > 0. The exponential term 
h(rij, rik) then provides the necessary asymptotics of the 
bond-order at large distances. In the above example, the 
contribution of bond i-k to fey will decay exponentially 
as rik increases. Usually, \x is treated as an adjustable 
parameter, but tight-binding bond theory tells us that 
for an expansion up to second moments and ignoring the 
contribution of 7r-orbitals the total energy needs to be 31 
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This is compatible with the empirical Tersoff-Brenner 




FIG. 3: (Color online) The bond-order 6y of bond i-j that 
sits at the surface of an exemplary diamond (110) surface 
depends on all neighbors. This includes neighbor k that sits 
on an opposite surface. At sufficient separation the influence 



of k on bij needs to vanish. 



formulation if we choose rj = X, 5 = 1/2 and 

2 



h(rij,r ik ) 



Pirn) 



(15) 



Using the functional form Eq. (4) for h(rij,ri k ) and 
Eq. (5) for j3(r) we obtain m = 1 and 



(16) 



Unfortunately, for m — I the value of /i contributes 
to the C44 shear modulus of the material. This is easily 
seen from the definition of this particular modulus: C44 
is given by 56 



a 



44 



1 d 2 E 
Wlh 2 



(17) 



where V is the volume of the crystal and E its total 
energy. The strain e characterizes the shear transforma- 
tion where all atoms are transformed from position r to 
f = (1 + T)r with 



e/2 e/2 
T= I e/2 e/2 
e/2 e/2 



(18) 



This particular transformation stretches some bonds in 
the diamond structure and contracts others. The second 
derivative of Eq. (17) then involves terms such as: 



d 2 h 



dri-jdnk 



m(m - l)(2fj li k) m (r i j - r ik ) m 2 h(r i j 1 r lk ) 



- m 2 (2 f i. lk ) 2m (r ZJ - r ik f m - A h{nj,r ik ) 



2m-2 1 



(19) 



In the equilibrium diamond structure = and this 
derivative vanishes only if m > 2. Choosing m — 1 would 
hence require a complete readjustment of all parameters 
to a set of material properties. For small deviations from 
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the crystalline ground-state, m = 3 hence removes the 
contribution of fi to the energy. While this is not fully 
consistent with Eq. (15), we use m = 3 in the following 
for convenience and to avoid refitting the potential. Since 
\i needs to have units of inverse length, we empirically 
choose [iik = Tq 1 to be the inverse of the dimer length ro 
of elements i-k. 

The Silicon potential of Kumagai et al. has a value 
of /i ^ that is independently fit. Here we therefor 

retain m = 1. Note that Kumagai et al. fit fi = 1.8 A 
while from Eq. (16) we obtain a value of fj, = 1.4 A 
We also note here that in our earlier screened REB02 
potential we enforced the proper limiting behavior for 
bij by an additional cutoff function /i(r\j, r^) — fci r ik) 
that depended on the distance nk only. For the potentials 
presented in this paper, we use the expression given by 
Eq. (4) because we believe that choosing a functional 
form close to that given by tight-binding bond theory is 
crucial for the transferability of the interatomic potential. 
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TABLE I: Parameters for the screened Tersolf (Till), Erhart 
& Albe (EA) and Kumagai potentials. 



C. Computational considerations 

A full cut-off free formulation as presented in the pre- 
ceding chapters is possible by computing a Voronoi tessel- 
lation of the atomic configuration in each time step. The 
screening functions would then be computed for atoms 
whose respective Voronoi cells share a face. However, 
this approach is computationally expensive and not linear 
scaling. In all practical cases, we therefor smoothly cut 
the interaction off at a certain distance to be able to use 
the usual linear scaling linked cell algorithms. 57 If this 
distance is large, the modulation of the bond-integrals 
will be weak and their asymptotic behavior essentially 
conserved. The final expression for the switching func- 
tion we use is hence 

Si, = fsim) + (1 - /s(r«))/c(r«)S« (20) 
with fc(r) = fs(r) that switches between radii r* and 

The specific parameters for the potentials presented in 
this article are given in Tab. I. The parameters are cho- 
sen with the following considerations in mind: n and r-i 
must lie between the first and second neighbor shell in 
the diamond or 3C structure (for C, Si and Si-C) and 
between the first and second neighbor shell in graphite 
(for C). Furthermore, r\ and r 2 for Si-Si must be smaller 
than the first Si-Si neighbor shell in 3C Si-C. The lat- 
ter constraint is the reason why r\ and r 2 for Si-Si are 
smaller than for C-C and Si-C if compared to the crys- 
talline bulk bond length. The outermost cutoff r% must 
be large enough to eliminate spurious peaks in the dimer 
force curves and the cohesive stress functions discussed 
below. This is usually achieved at about r% « 2.5r nn 
where r nn is the nearest neighbor distance in the dia- 
mond or 3C structure. We furthermore empirically fix 
r 2 = 1.2n and r\ — 2r*. For the Tersoff potential we 



use the original Tcrsoff-Lorentz-Berthclot 4 mixing rule 
^SiC = V rcrsi for n and r*. The values of C m i n and 
C max are chosen such that for three atoms located on the 
corners of an equilateral triangle three unscreened bonds 
exist, and for four atoms on the corners of a square four 
bonds exist. 



IV. PROPERTIES OF THE SCREENED 
POTENTIALS 

We report some select properties of the screened po- 
tentials and compare those to their unscreened counter- 
parts and higher level quantum calculations. 58 In what 
follows, we denote Tersoff's third-generation potential 4 
as Till, and the screened incarnation as TIII+S. Simi- 
larly, we denote Erhart & Albe's potential 13 as EA, and 
the screened incarnation as EA+S. Kumagai et al.'s 49 po- 
tential will be referred to as Kumagai and Kumagai+S 
in it's unscreened and screened incarnation, respectively. 
For completeness, we also compare to results obtained 
with the REB02 10 and screened REB02 (REB02+S) 36 
potential for carbon, and the Stillinger- Weber (SW) po- 
tential for silicon. 

If not otherwise noted, DFT reference calculations are 
carried out by us and employ the local density approx- 
imation 61 and projector augmented waves. 62 The wave 
functions are expanded on a real space grid. We use the 
GPAW code. 63,64 Table II lists some properties of dia- 
mond, silicon and 3C silicon carbide as obtained from 
the classical potentials and this particular DFT method. 
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240 


C2i (GPa) 






273" 


279 


311 


305 


7{m} (J m " 








4.17 9 


1.85 9 


1.67 9 


7{no} (J m" 


- 2 ) 






3.29 


2.40 


2.29 


7{ioo} (J m" 


- 2 ) 






5.46 9 


4.12 (4.21) 9 


3.87 (3.93) 9 


7 { 2 ioo} ( J m " 








3.48 9 


2.87 (2.85) 9 


2.92 (2.85) 9 



TABLE II: Properties of diamond, silicon and 3C silicon-carbide. Unless referenced, DFT results are LDA (see text). Values 
in parenthesis are for the unscreened potential, if different from their screened counterpart. C44 is the C44 modulus obtained 
without relaxation of atomic positions. The (111) surface is cut at the shuffle plane. a Ref. 65 6 Ref. 66 c Ref. 67 d Fracture energy, 
Ref. 68 e Ref. 69 ^DFT-GGA, Ref. 70 9 Ref. 71 h Refs. 72-74 *Ref. 75 j Ref. 76 fe Ref. 77 'Ref. 78 m Ref. 79 "Ref. 80 °The Si (110) 
surface reconstructs in DFT-LDA. The empirical potentials do not capture this reconstruction. p Fixed to the experimental 
value. 'Energies obtained by creating a silicon and a carbon terminated surface. 
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A. Fracture 

1. Cohesive stress 

We compute the cohesive stress functions by separat- 
ing the ideal bulk of the crystal for diamond, silicon, and 
3C silicon-carbide to create (100), (110) and (111) sur- 
faces. These calculations are carried out unrelaxed, and 
the cohesive stress function that is shown in Fig. 4 is 
the first derivative of the total energy curves obtained 
normalized by the exposed surface area. The value of x 
denotes the distance of the newly created surface such 
that x = is the limit of the bulk crystal and x — > oo 
are two free surfaces. 

The maximum force obtained for all structures and all 
surfaces probed here is in reasonable agreement with the 
DFT calculations. However, the asymptotic behavior of 
the cohesive stress is significantly lower than the values 
obtained from DFT calculations for the (110) and even 
worse for the (111) surface. Since the surface energy is 
the area beneath the cohesive stress functions of Fig. 4, 
this difference can be attributed solely to a mismatch in 
surface energy. We list the energies for these surfaces in 
Tab. II. While all potentials give reasonable values for 
the high energy (100) surface, the agreement with DFT 
calculation for the (110) and (111) surfaces are worse. 
For carbon and 3C silicon-carbide, the order of the en- 
ergetics of (110) and (111) surfaces is reversed in DFT 
calculations and experiments. All potential except for 
the REB02+S follow the experimental order. It is some- 
what surprising that these classical potentials appear to 
capture the peak force at the transition state more accu- 
rately than the equilibrium surface energies. 

For opening a (100) diamond surface we find that the 
force for TIII+S and EA+S has two distinct peaks, with 
the peak at the larger separation having a higher force. 
These peaks are less pronounced, but also visible, for 
the silicon- carbide (100) surface, but do not show up on 
the (110) or (111) surfaces. The origin of this is a too 
sudden drop in Xij that stems from choosing m = 3, and 
not m = 1 in Eq. (3) as the rigorous bond-order theory 
suggests. 31 This in return leads to an overestimation of 
bij for the transition state and hence a potential that is 
too attractive in that region. 

Finally, we note that the potential energy landscape of 
the REB02+S is more corrugated than the one obtained 
for TIII+S, EA+S and Kumagai+S. This is related to the 
treatment of 7r-electron in the REBO formalism. In brief, 
an additive correction is applied to bij and Xij (given in 
Eqs. (2) and (3), respectively). The value of that correc- 
tion depends on the coordination numbers of the atoms 
in the vicinity of the bond and was fit to the atomization 
energies of a select set of hydrocarbon molecules. Since 
coordination numbers are integer values, the transition 
values upon changes in coordination are obtained from a 
cubic spline interpolation. This cubic spline is the origin 
of the additional corrugation seen for the REB02+S in 
Fig. 4. We also note that for (110) and (111) surfaces 



the coordination number jumps from 4 for a bulk atom 
to 3 for a surface atom. On the (100) surface, the co- 
ordination number jumps from 4 in the bulk to 2 at the 
surface giving rise to an additional transition state with 
coordination number 3 that is the origin of the peaks 
seen in Fig. 4 for REB02+S on this particular surface. 
The simpler formulation given by Eqs. (1) to (8) without 
the spline corrections that is the basis of the Till, EA 
and Kumagai potentials has the advantage that it yields 
a smoother potential energy landscape, albeit at the cost 
of limited accuracy in particular in complex molecular 
systems. 

2. Static crack 

In addition to the cohesive stress function we com- 
pute bond-breaking events in a mode I crack geometry 
using the method by Perez and Gumbsch. 37,38 In brief, 
we consider a small atomistic region around the crack 
tip and fix the boundary atoms of this region using the 
near field solution of the displacements from linear elas- 
tic fracture mechanics. Then, the stress intensity factor 
K is increased step-wise, the system is relaxed, 81 and we 
monitor the length of the bond in front of the crack tip. 
We also investigate the closing of a crack by decreasing 
the stress intensity factor and monitoring the length of 
the bond behind the crack tip. In all calculations the 
crack tip is centered on the bond of interest. More in- 
formation on the technique can be found in Refs. 37 and 
38. 

Results for a crack on the (110) surface with a [110] 
crack front for diamond, silicon and 3C silicon-carbide are 
shown in Fig. 5. We do not show the unscreened poten- 
tials which do not break bonds in this kind of simulation. 
For TIII+S and EA+S the agreement with DFT calcula- 
tions is reasonable. For diamond and silicon, the TIII+S 
follows the DFT results almost exactly in predicting the 
correct stress intensity factor K for bond breaking and 
bond formation. EA+S overestimates the stress inten- 
sity factor K + required for breaking and underestimates 
the stress intensity K_ for bond formation hence giving 
a too large lattice trapping region AK = K + — K_ for 
carbon. For silicon, the width of the lattice trapping re- 
gion AK is well described by both potentials. For 3C 
silicon-carbide, TIII+S and EA+S give almost identical 
results but overestimate the lattice trapping AK. Addi- 
tionally, the opening of the bond in our DFT calculations 
proceeds more smoothly. This could be related to charge 
transfer that occurs in silicon-carbide and is not captured 
by our potentials. 

B. Melting 

We determine the melting point for diamond, silicon 
and 3C silicon-carbide by equilibrating a crystal-melt in- 
terface in a simulation without heat exchange with some 
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FIG. 4: (Color online) Cohesive stress a (derivative of the energies obtained by separating a bulk crystal normal to certain 
surfaces normalized by their area) for carbon, silicon and 3C silicon-carbide. Here we show this function for the creation of 
low-index (100), (110) and (111) surfaces. Broken lines show results for the unscreened interatomic potentials. Solid lines show 
their screened counterparts and DFT results. Curves are shifted vertically to be distinguishable. The (111) surface is cut at 
the shuffle plane. 



external bath. In these simulations, the crystal-melt in- 
terface advances or recedes until the system is equili- 
brated to the melting temperature. In all cases the (100) 
surface is exposed to the melt and the pressure is con- 
trolled by an anisotropic Andersen barostat 82 that con- 
trols the box size independently in all three Cartesian 
directions. 

The result of this calculation are summarized in 
Tab. III. The melting points for diamond are taken at 
the pressure of the diamond/graphite/melt triple point 



(12.4 GPa) and lie in the range of experimental values 
for all screened potentials. We note that while 4713 K 
seems to be the universally referenced melting point of 
diamond, the experimental values spread over a much 
larger range with initial melting reported at temperatures 
as low as 3500 K. 83 No melting point could be obtained 
for the unscreened potential because bulk diamond spon- 
taneously transforms into a graphite under these pres- 
sure/temperature conditions. Since the interaction range 
of the unscreened case is considerably smaller than the 



9 



2.5 



2.0 



1.5 



i n 



Tersoff+S 
Erhart & Albe+S 
Kumagai+S 
DFT-LDA 




C" 



Si 



3C-Si-C 



Expt. 

tiii 

TIII+S 

EA 
EA+S 
REB02 
REB02+S 



3500— 5000 6 , 4713 c 1687= 2818 ± 40 e 



5240 ± 30 

_d 

4210 ± 25 

_d 

3950 ± 20 



2580 ± 15 

2330 ± 15 3190 ± 15 
2510 ± 15 

2365 ± 15 3235 ± 15 




0.6 0.8 1.0 1.2 1.4 



K/K G 

FIG. 5: (Color online) Bond length r before the crack tip (left 
three curves) and after the crack tip (right three curves) for a 
crack on the (110) surface with a [110] crack front in carbon, 
silicon and 3C silicon-carbide as a function of the stress inten- 
sity factor K. The bond length r is here scaled by the bond 
length ro in the bulk crystal, and the stress intensity factor 
is scaled by the factor Kg obtained from Griffith's criterion 
using relaxed surface energies. 



interlayer graphite spacing this conversion can proceed 
without a volume expansion and hence without perform- 
ing work against the external pressure. The screened 
potential have a longer range. The individual graphitic 
sheets do interact and inhibit this transition at suffi- 
ciently high pressures. 

The melting point for silicon at zero pressure is over- 
estimated by about 1000 K by both the TIII and EA po- 
tentials. This overestimation has been noted before, 13,85 
and Kumagai and co-workers pointed out that it is re- 
lated to the angular term. The Kumagai potential em- 
ploys a different angular term and does correct the melt- 
ing point as independently confirmed by Schelling 86 and 
here. However, the improved melting point comes at an 
expense of surface energies that are considerably lower 



Kumagai - 1725 ± 10 

Kumagai+S - 1625 ± 10 

SW - 1636 ±5 f 

TABLE III: Melting points of crystalline diamond, silicon 3C- 
Si-C in Kelvin for the different potentials studied here. The 
melting point was determined by equilibrating a (100) surface 
with the melt. Error is the standard deviation of the temper- 
ature fluctuation in an equilibrated NVE ensemble of 18432 
atoms. a At 12.4 GPa. 6 Ref. 83 c Ref. 66 d A diamond/melt 
interface is unstable in the unscreened potentials, see text. 
e Ref. 84 f System size of 65536 atoms total. 



than Till and EA energies which themselves are an un- 
derestimation of the respective DFT results (see Tab. II). 
For silicon-carbide we obtain melting points that are only 
about 400 K too high. TIII and EA solids melt at roughly 
identical temperatures. In all cases, the screening func- 
tion lowers the melting point compared to the respective 
unscreened potential by about 10 to 15%. 



C. Glass formation 

1. Hybridization of amorphous carbon 

Classical empirical bond-order potentials notoriously 
fail at describing the properties of amorphous carbon 
that is quenched from the melt. One particular property 
that is also accessible from experiments and ab-initio cal- 
culations it the fraction of diamond-like, sp 3 hybridized 
atoms as a function of the density of the amorphous sam- 
ple. For example, the Tersoff, REBO, 5,6 and REB02 po- 
tentials are known to fail to describe this relationship and 
typically yield 20% to 40% of sp 3 close the the density of 
diamond where the sp 3 fraction should saturate. 36,87 ' 88 
For deposition processes, a common cure is to slightly 
increase the cut-off range of the potential but keeping 
it between the first and second nearest neighbor shell of 
graphite and diamond. 89,90 This cure only works above 
a certain density. 14 

Here, we compute sp 3 (p) curves by quenching liquid 
carbon within 0.5 ps from 5000 K to 300 K at constant 
volume. The same procedure has been used in ab-initio 91 
and non-orthogonal tight-binding (NOTB) 36 calculations 
that will be used as a reference here. We also report the 
experimental analysis of physically deposited amorphous 
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FIG. 6: (Color online) Carbon: Fraction of sp 3 in amorphous 
carbon that is quenched from a 5000 K melt to 300 K with a 
time constant of 0.5 ps. Experimental data is from Ref. 92, 
DFT data is from Ref. 91 and NOTB data is from Ref. 36. 



carbon of Ref. 92 for comparison. An atom contributes 
towards the sp 3 fraction if it has four neighbors within a 
distance of 1.85 A. 

All this data, along with results for the screened and 
unscreened potentials discussed in this work are shown in 
Fig. 6. The TTII+S potential follows the NOTB data al- 
most exactly, albeit yielding an sp 3 fraction that is lower 
by a few percent. EA+S also follow the NOTB curve, 
but the sp 3 fraction is lower than the one obtained by 
TIII+S. All unscreened potentials are worse, predicting 
at best 50% to 60% sp 3 at densities of 3.5 g cm -3 where 
the sp 3 fraction should be around 80%. 
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FIG. 7: (Color online) Silicon: Volume per atom as a function 
of temperature when quenching from the liquid phase to 300 K 
at a rate of 1 Kps -1 and zero external pressure. The periodic 
cell contained 4001 atoms in all cases. Experimental data is 
taken from Ref. 97. 



liquid phase density consistent with experiments. 97 The 
screened and unscreened Kumagai potential give simi- 
lar results. The screened Kumagai however seems to be 
marginally better at reproducing the slope dV/dT of the 
experimental temperature dependence that was reported 
by Rhim et el. 97 

The temperature at which the density peaks during 
solidification is typically associated with the glass transi- 
tion temperature T g . 95 Both Kumagai and SW potentials 
give a T g of about 1000 K in excellent agreement with 
measurements. 98 Till and EA overestimate both glass 
transition and melting temperature by about a factor of 
1.5. 



3. Pair distribution functions of amorphous silicon- carbide 



2. Supercooling amorphous silicon from the melt 

In computer simulations, amorphous silicon is typically 
quenched from the melt at constant pressure rather than 
constant volume. 93-95 We here carry out such simulation 
at zero external pressure and quench rates of 1 Kps -1 us- 
ing Berendsen temperature and pressure control 96 with 
relaxation time constants of approximately 1 ps for tem- 
perature and 10 ps for pressure. The quench starts from 
the melt equilibrated at 3000 K. 

We first note that the density of the melt does no- 
tably depend on the potential under consideration. Fig. 7 
shows the atomic volume V as a function of temperature 
T during the quench. The volume at the highest tem- 
perature (3000 K) is the equilibrated melt. All potentials 
but the screened TIII+S predict a liquid phase that is 
denser than the supercooled amorphous that is shown at 
300 K. However, only Kumagai, Kumagai+S and the SW 
potential predict a liquid phase that is denser than the 
crystalline. The densest liquid phase is given by the Ku- 
magai potential which is the only potential to reproduce a 



Finally, we also report pair distribution functions of 
quenched amorphous silicon carbide. Silicon-carbide is 
quenched at zero external pressure using the procedure 
outlined in the previous section for silicon. Figure 8 
summarizes the results alongside experimental data from 
Ref. 99. All potentials reproduce the experimental pair 
distribution functions reasonably. The unscreened po- 
tentials give pair distribution functions that are essen- 
tially indistinguishable from their screened counterparts 
and therefor not shown. In all cases, the experimental 
data is broader than the data obtained from our simu- 
lations. This is probably attributable to additional line 
broadening mechanisms that are active in the respective 
experimental setup. Also, the experimental amorphous 
Si-C was created by ion irradiation and not by quench- 
ing, which could be the origin of some of the observed 
differences. 

The notable differences between the two potentials are 
the heights of the nearest-neighbor peaks. The TIII+S 
potential overestimates the height of the first neighbor 
peak significantly. From the distribution functions for 
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FIG. 8: (Color online) Silicon-carbide: Pair distribution func- 
tions for stochiometric silicon carbide. The final amorphous 
structure was obtained by quenching a box of 4000 atoms 
from the melt to 300 K at a rate of 1 K ps _1 and zero external 
pressure. Curves are shifted vertically to be distinguishable. 
The experimental data is taken from Ref. 99. 



pure amorphous carbon (not shown) we see that this 
peak corresponds to the C-C bond length. The TIII+S 
potential also appears to overestimate the peak at 2.5 A 
that is barely visible in the EA+S simulation and the 
experimental data. This length corresponds roughly to 
the Si-Si bond lengths. Hence, the TIII+S appears to 
favor dimerization over the formation of a homogeneous 
melt, leading to a somewhat different structure than that 
found in experiments. 

V. CONCLUSIONS 

We have presented a simple method to augment exist- 
ing bond-order potential by changing their cut-off pro- 



cedure. This fixes a number of issues with the descrip- 
tion of non-equilibrium properties of matter, such as frac- 
ture or amorphous phase formation. We here stress that 
without any reparameterization of the potentials we are 
able to obtain correct cohesive stresses, proper bond- 
breaking in mode I cracks and appropriate properties of 
the amorphous phase. Both the Tersoff III and Erhart 
& Albe's potential are fitted to ground-state properties, 
yet they are able to reasonably describe these transition 
states. The potential energy expression given by Eq. (1) 
is hence an exquisite extrapolation scheme. Surely, this 
is due to the fact that there are good theoretical argu- 
ments 2,8 ' 31 ' 32,100 for this particular functional form. Fu- 
ture work will focus on augmenting a recent potential for 
the ternary Si-C-H system in a similar manner. 19 



Force routines for the potentials of this paper are avail- 
able at the location given in Ref. 58. 
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